/*  sdfcoshs.cpp  by K.Tsuru  */
// function ID 3307     DRADIX
#ifndef SN_H
#include "sn.h"
#endif
/**************************************************************
SDouble class
hyperbolic function cosh x by series for |x| < DRADIX
If first != 0 calculate cosh(x), else cosh(x) - 1 for small x.
Default argument is declared as first = 1.
When x is larger than about 20.0, it is slow and error is large.
***************************************************************/
static const char* func = "CoshSeries";
SDouble CoshSeries(const SDouble& x){
    return CoshSeries(x, 1);
}
SDouble CoshSeries(const SDouble& x, int first){
    if(x.NetRdxExp() > 1) x.SetError(x.OUT_OF_RANGE, func, 3307);

    SDouble sum, xsq = x*x, delta, t;
    if(first){
        sum = 1.0;          sum.FixedPoint(sum.RdxExp());
    }else{
        sum.SetZero();      sum.FixedPoint(xsq.RdxExp());
    }
    delta = DsDiv(xsq, 2uL);    // delta = x*x/2

    ulong n = 3uL, den, mt = sum.SlOpMaxValue();
    sum = DDAdd(sum, delta);    // sum += delta;
    while( delta.Sign(3307) ){
        den = n*(n+1uL);
        if(den >= mt) break;
        t = delta*xsq;
        delta = DsDiv(t, den);
        sum = DDAdd(sum, delta);
        n += 2uL;
    }

    while( delta.Sign(3307) ){
        if(n >= mt){
            sum.SetError(sum.NOT_CONVERGE, func, -3307);
            break;
        }
        t = delta*xsq;
        delta = DsDiv(t, n);    delta = DsDiv(delta, n+1uL);
        sum = DDAdd(sum, delta);
        n += 2uL;
    }
    sum.PointFree();
    sum.Reform(3307);
    sum.upToTerm = (n+1)/2;
    return sum;
}